function L = graphLaplacian(S,varargin)
    par.normalize = 'unnorm';
    par = process_parameter(par, varargin{:});
    
    [m,n] = size(S);

    switch par.normalize
      case 'unnorm'   
        D = spdiags(sum(S,2),0,m,n);
        L = D - S;
      case 'symnorm'
        D = spdiags(sum(S,2).^(-0.5),0,m,n);
        L = spdiags(ones(m,1),0,m,n) - D * S * D;
      case 'rwnorm'
        D = spdiags(1./sum(S,2),0,m,n);
        L = spdiags(ones(m,1),0,m,n) - D*S;
      otherwise
        error('unknown normalization type')
    end
end
